Feedback-optimized parallel tempering Monte Carlo 
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We introduce an algorithm to systematically improve the efficiency of parallel tempering Monte 
Carlo simulations by optimizing the simulated temperature set. Our approach is closely related to a 
recently introduced adaptive algorithm that optimizes the simulated statistical ensemble in general- 
ized broad-histogram Monte Carlo simulations. Conventionally, a temperature set is chosen in such 
a way that the acceptance rates for replica swaps between adjacent temperatures are independent 
of the temperature and large enough to ensure frequent swaps. In this paper, we show that by 
choosing the temperatures with a modified version of the optimized ensemble feedback method we 
can minimize the round-trip times between the lowest and highest temperatures which effectively in- 
creases the efficiency of the parallel tempering algorithm. In particular, the density of temperatures 
in the optimized temperature set increases at the "bottlenecks" of the simulation, such as phase 
transitions. In turn, the acceptance rates are now temperature dependent in the optimized tem- 
perature ensemble. We illustrate the feedback-optimized parallel tempering algorithm by studying 
the two-dimensional Ising ferromagnet and the two-dimensional fully-frustrated Ising model, and 
brieffy discuss possible feedback schemes for systems that require configurational averages, such as 
spin glasses. 

PACS numbers: 75.50.Lk, 75.40.Mg, 05.50.+q 



I. INTRODUCTION 



The free energy landscapes of complex systems are 
characterized by many local minima that are separated 
by entropic barriers. The simulation of such systems with 
conventional Monte Carlo methods is slowed down by 
long relaxation times due to the suppressed tunneling 
through these barriers. Extended ensemble simulations 
address this problem by broadening the range of phase 
space which is sampled in the respective reaction coordi- 
nate. Recently, an adaptive algorithm I'l has been intro- 
duced that explores entropic barriers by sampling a broad 
histogram in a reaction coordinate and iteratively opti- 
mizes the simulated statistical ensemble defined in the re- 
action coordinate to speed up equilibration. The key idea 
of the approach is to measure the local diffusivity along 
the reaction coordinate, thereby identifying the bottle- 
necks in the simulations and then using this information 
to systematically shift statistical weights towards these 
bottlenecks in a feedback loop. The optimized histogram 
converges to a characteristic form exhibiting peaks at the 
bottlenecks of the simulation, e.g., in the vicinity of the 
entropic barriers. The simulation of an optimized ensem- 
ble results in equilibration times that can be substantially 
lower compared to other extended ensemble simulations 
that aim at sampling a flat histogram in the respective 
reaction coordinate Flat-histogram algorithms in- 

clude the multicanonical method '5, 'ij , broad histograms 
[3 and transition matrix Monte Carlo 6] when combined 
with entropic samplingj_as_well as the adaptive algorithm 
of Wang and Landau 



Parallel tempering Monte Carlo has proven to be 
a strong "workhorse" in fields as diverse as chemistry, 
physics, biology, engineering, and material science [lOl |. 
Similar to replica Monte Carlo Tl'l , simulated tempering 
[T^ . or extended ensemble methods JJI], the algorithm 
aims to overcome entropic barriers in the free energy 
landscape by simulating a broad range of temperatures. 
This allows the system to escape metastable states when 
wandering to higher temperatures and subsequently re- 
laxing at lower temperatures again on time scales that 
are substantially smaller than conventional simulations 
at a fixed temperature. In this paper, we maximize the 
efficiency of parallel tempering Monte Carlo by optimiz- 
ing the distribution of temperature points in the simu- 
lated temperature set such that round-trip rates of repli- 
cas between the two extremal temperatures in the sim- 
ulated temperature set are maximized. The optimized 
temperature sets are determined by an iterative feed- 
back algorithm that is closely related to the previously 
mentioned adaptive ensemble optimization method for 
broad- histogram Monte Carlo simulations jljj . The feed- 
back method concentrates temperature points near the 
bottleneck of a simulation, typically in the vicinity of a 
phase transition or the ground state of the system. As a 
consequence, we find that for the optimal choice of tem- 
peratures the acceptance probabilities for swap moves 
between neighboring temperature points show a strong 
modulation with temperature and are not independent 
of temperature as suggested in several recent approaches 

niiaiiiiiiiiiia. 

The paper is structured as follows: In Sec^Jwe present 
a detailed introduction of the parallel tempering Monte 



Carlo method. In Sec. IIIII we generalize the feedback 
method of Ref. |1| to the parallel tempering Monte Carlo 
algorithm. Results on two paradigmatic models, the two- 
dimensional Ising ferromagnet and the two-dimensional 
fully-frustrated Ising model are presented in Sec. IIVI as 
well as a discussion on how to proceed with systems that 
require configurational averages, such as spin glasses. 
Concluding remarks follow in Sec. 

II. PARALLEL TEMPERING MONTE CARLO 

In the parallel tempering Monte Carlo algorithm ^2;, 
ITU IT^ IT^ . M non- interacting replicas of the system 
are simulated simultaneously at a range of temperatures 
{Ti, T2, . . . , Tm}- After a fixed number of Monte Carlo 
sweeps a sequence of swap moves, the exchange of two 
replicas at neighboring temperatures, Tt and Ti+i, is sug- 
gested and accepted with a probability 

p(£;„r, T,+i)=min{l,exp(A/3Ai;)} , (1) 

where A/3 = — l/T^ is the difference between the 

inverse temperatures and Ai? = Ei+i — Ei is the differ- 
ence in energy of the two replicas. At a given tempera- 
ture, an accepted swap move effects in a global update 
as the current configuration of the system is exchanged 
with a replica from a nearby temperature. For a given 
replica, the swap moves induce a random walk in tem- 
perature space. This random walk allows the replica to 
overcome free energy barriers by wandering to high tem- 
peratures where equilibration is rapid and return to low 
temperatures where relaxation times can be long. The 
simulated system can thereby efficiently explore complex 
energy landscapes that can be found in frustrated spin 
systems 20], spin glasses |2^ |23j or proteins [23 |. 
While the simulation of M replicas takes M times more 
CPU time, the speedup attained with parallel temper- 
ing Monte Carlo can be orders of magnitude larger. In 
addition, one often wishes to measure observables as a 
function of temperature. Thus parallel tempering Monte 
Carlo delivers already several measurements at different 
temperatures in one simulation. In order to maximize the 
number of statistically independent visits at low temper- 
atures, we want to maximize for each replica the number 
of round-trips between the lowest and highest tempera- 
ture, Ti and Tm, respectively. The rate of round-trips 
of a replica strongly depends on the simulated statisti- 
cal ensemble, that is the choice of temperature points 
{Ti, T2, . . . , Tm} in the parallel tempering simulation. 

In this paper, we present an algorithm that systemati- 
cally optimizes the simulated temperature set to max- 
imize the number of round-trips between the two ex- 
tremal temperatures Ti and Tf,[ for each replica and 
thereby substantially improve equilibration of the system 
at all temperatures. Conventional approaches assume 
that to achieve this goal, the simulated temperature set 
should be chosen in such a way that the probability for 
replica swap moves at neighboring temperatures should 
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FIG. 1: (Color online) Sketch of the random walk that a 
given replica performs in temperature space in the course of 
the simulation. Ideally, the replica will wander up (rup) and 
down (rdown) in the simulated temperature range [Tmin, Tmax]. 
The goal of the feedback method is to maximize the number of 
round trips each replica performs in this temperature range, 
and thereby minimize the average round-trip time Trt = Tup + 

^down . 

be "flat" , that is approximately independent of temper- 
ature. If the specific heat of the system is assumed to be 
constant, then a good approximation for such a temper- 
ature set can be attained with a geometric progression 
|17| . Given a temperature range [Ti , Tm\ , the intermedi- 
ate M — 2 temperatures can be computed via 



fc-i 

Tk=Til[ R, 




The geometric progression peaks the number of tem- 
peratures around the minimum temperature Ti where 
a slower relaxation is assumed. This is not optimal when 
the system has a diverging specific heat (at an interme- 
diate temperature): Because in order to ensure enough 
overlap between the energy distributions of neighboring 
temperatures ATi^^+i ^ CyTi/^/N, where Cy is the spe- 
cific heat per spin and N the number of spins, the accep- 
tance probabilities are inversely correlated to the func- 
tional behavior of Cy via the inverse beta function law 
|l7j . Thus, for example at a phase transition where the 
specific heat diverges, the acceptance probabilities for a 
geometric temperature set will show a pronounced dip 
(cf. Sec. II V A|) . More complex methods such as the ap- 
proach by Kofke [l^ its improvement by Rathore 
et a l. .16.1 . as well as the method suggested by Predescu 
[17L llSjl aim to obtain acceptance probabilities for the 
parallel tempering moves that are independent of tem- 
perature by compensating for the effects of the specific 
heat. In particular, Kone and Kofke suggest that an 
acceptance probability of 23% is optimal. In this work 
we show that this is not necessarily the optimal case. 



III. TEMPERATURE SET OPTIMIZATION 

Our goal is to vary the temperature set {Ti} of a par- 
allel tempering simulation in such a way that for a given 
system we speed up equilibration at all temperatures. To 
accomplish this, we maximize the rate of round trips that 
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each replica performs between the two extremal temper- 
atures Tmin = Ti and Tmax = Tm following a similar ap- 
proach to the ensemble optimization technique presented 
in Ref. I*]. For a given temperature set, we can measure 
the diffusion of a replica through temperature space by 
adding a label "up" or "down" to the replica that indi- 
cates which of the two extremal temperatures, Tmin or 
Tmax respectively, the replica has visited most recently. 
The label of a replica changes only when the replica visits 
the opposite extreme. For instance, the label of a replica 
with label "up" remains unchanged if the replica returns 
to the lowest temperature Tmin, but changes to "down" 
upon its first visit to Tmax- This is illustrated in Fig. ^ 
For each temperature point in the temperature set {Ti}, 
we record two histograms 7iup(Ti) and ndown{Ti). Be- 
fore attempting a sequence of swap moves, we increment 
that histogram at temperature Ti which has the label of 
the respective replica currently at temperature Ti. If a 
replica has not yet visited either of the two extremal tem- 
peratures, we increment neither of the histograms. This 
allows us to evaluate for each temperature point the frac- 
tion of replicas which have visited one of the two extremal 
temperatures most recently (e.g., T„iin) as 



riupiTi) 



^down 



(3) 



The labeled replicas define a steady-state current j 
from Tiiiiii to Tmax that is independent of temperature, 
e.g., the rate at which "up" walkers are created at Tmin 
and - in equilibrium - absorbed at Tmax- In the following 
we assume that T is a continuous variable, independent 
of the temperature points in the current temperature set. 
We can then determine the current j to first order in the 
derivative as 



J = D{T)rjiT) 



dT 



(4) 



where D{T) is the local diffusivity at temperature T and 
the derivative df /dT is estimated by a linear regression 
based on the measurements in Eq. il{T) is a density 
distribution indicating the probability for a replica to 
reside at temperature T. We approximate ri{T) with a 
step-function ri{T) = C/AT, where AT T^+i-T, is the 
length of the temperature interval around temperature 
Ti < T < Ti+i for the current temperature set. The 
normalization constant C is chosen such that 



7]iT)dT = C / — = 1 . 
T, ^ ' Jt, at 



(5) 



Rearranging Eq. gives a simple measure of the local 
diffusivity D{T) of a replica at temperature T 



D{T) 



AT 
dfJdT 



(6) 



where we have dropped the normalization constant C and 
the current j which is constant for any specific choice of 
temperature set. 



To increase the efficiency of the algorithm, we maxi- 
mize the current j in temperature space by varying the 
simulated temperature set {Ti} and thus varying the 
probability distribution r]{T) between the two extremal 
temperatures, Tmin and Tmax, which are not changed. In 
Ref. 1 it has been shown that the optimal probability dis- 
tribution r/°P'(T) is inversely proportional to the square 
root of the local diffusivity D{T): 



r/°P*(T) cx 



1 



(7) 



For the optimal distribution of temperature points the 
fraction /°P'(r) then decays as 



df °P' 1 
dT J « ^yopt ' 



(8) 



which implies that for any given temperature interval 
AT = Tijf.1 — Ti of the optimal temperature set the frac- 
tion has a constant decay 

A/opt = jopt(y^) _ /opt(j.^^) ^ ^/(^^ _ ^) ^ 

where M is the number of replicas. 

In our algorithm we approach the optimal temperature 
set and its respective probability distribution by itcra- 
tively feeding back measurements of the local diffusivity. 
After measuring the diffusion of replicas for a given tem- 
perature set an improved probability distribution rj'iT) 
is found as 



c 

AT' 



C" 



J_dl 
AT dT 



(10) 



where the normalization constant C is again chosen so 
that the normalization condition in Eq. jSJl is met. The 
step- function r]'(T) is still defined for the ori^ma/ temper- 
ature set, that is the jumps in the function occur at the 
temperature points in {Ti}. The optimized temperature 
set {r/} is then found by choosing the k-th temperature 
point such that 



Ti 



T]'{T)dT 



M ' 



(11) 



where 1 < k < M and the two extremal temperatures 
Tl — Ti and Tjjj = Tm remain fixed. 

We summarize the feedback algorithm by the following 
sequence of steps 

• Start with a trial temperature set {Ti}. 

• Repeat 

o Reset the histograms nup(T) = ndown(T) — 0. 

o For the current temperature set perform a 
parallel tempering simulation with Ngyj swap 
moves. After each sequence of swap moves: 

Update the labels of all replicas. 
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Record histograms nup(T) and ndown(T'). 

o For the given temperature set estimate an op- 
timized probabiUty distribution rj' {T) via 

o Obtain the optimized temperatures {T[} via 

o Increase the number of swaps N^vi *— 2A^sw 

• Stop once the temperature set {T^} has converged. 

The initial number of swaps Ng^ should be chosen large 
enough such that a few of round-trips are recorded, 
thereby ensuring that steady-state data for nup(T) and 
JT-downlr) are measured. The derivative df/dT can be 
determined by a linear regression, where the number 
of regression points is flexible. Initial batches with the 
limited statistics of only a few round trips may require 
a larger number of regression points than subsequent 
batches with smaller round-trip times and better statis- 
tics. 



IV. RESULTS 

A. Ferromagnetic Ising model 

In order to illustrate the feedback method, we start 
with a simple test model, the two-dimensional ferromag- 
netic Ising model (FM). The Hamiltonian for the model 
is given by 

Hfm — — J ^ SiSj , (12) 

where J = 1 and Si = ±1 represent Ising spins on 
a square lattice with N = L"^ spins. In our simula- 
tions we apply periodic boundary conditions and con- 
sider nearest-neighbor interactions only. The simple Ising 
model with uniform couplings has no frustration or dis- 
order, and exhibits a second-order phase transition at 
Tc = 2/ln(l V2) ~ 2.269 from a magnetically ordered 
to a paramagnetic phase. 

For simplicity, we define an initial temperature set 
{Ti\ with M = 21 temperature points performing a ge- 
ometric progression, Eq. for a temperature interval 
\Ti — 0.1, T^v/ = 10.0]. The minimum temperature Ti 
is chosen low enough such that the system can approach 
the zero-temperature ground state of the model, and the 
maximum temperature Tm is chosen well above the crit- 
ical region of the phase transition. For a short paral- 
lel tempering simulation (A^sw = 1.6 • 10^, one parallel 
tempering swap after each lattice sweep) using this ini- 
tial temperature set, we measure the diffusive current of 




T [index] 

FIG. 2: (Color online) Fraction /(T) of replicas diffusing from 
the lowest to the highest temperature as a function of the 
temperature index for the ferromagnetic Ising model. For 
the initial temperature set based on a geometric progression 
(filled squares), the fraction shows a sharp drop between two 
temperature points. A similar behavior is found for a tem- 
perature set where the acceptance rates are constant ~ 40% 
independent of temperature (temperature set with "flat" ac- 
ceptance rates, open squares). In contrast, for the optimized 
temperature set (triangles) the fraction constantly decreases. 
The inset shows the fraction as a function of temperature T. 
The dashed line in the inset represents the critical tempera- 
ture of the two-dimensional Ising model, Tc f» 2.269. 



replicas wandering from the lowest to the highest tem- 
perature using an additional label for each replica as de- 
scribed above. In Fig.|2we show the fraction of replicas 
diffusing "up" in temperature space. For the geomet- 
ric temperature set a sharp drop between two tempera- 
ture points is observed close to the critical region of the 
phase transition at Tc « 2.269. Similar results are also 
found for a "flat" temperature set where the acceptance 
rates are approximately constant around ^ 40% (with 
fluctuations of ~ 10%) and independent of temperature, 
although the drop is not as pronounced. 

Calculating the local diffusivity D{T) for the random 
walk in temperature space using Eq. (|SJ| , we find a strong 
suppression around this critical region as illustrated in 
Fig. 121 When increasing the size of the simulated sys- 
tem, the dip in the local diffusivity further proliferates, 
an additional indicator that the slowdown of the random 
walk in temperature space is dominated by the occur- 
rence of a phase transition. Note the logarithmic scale of 
the ordinate axis in Fig. O 

When applying the feedback, Eq. lll|l . to define a new 
temperature set, this suppression in the local diffusivity 
leads to a concentration of temperature points near the 
critical temperature as shown in Fig. ^ The feedback 
thereby tries to compensate for the reduced mobility of 
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FIG. 3: (Color online) Local diffusivity D{T) of the random 
walk a replica performs in temperature space for the ferro- 
magnetic Ising model as a function of temperature T after 
the feedback optimization for several system sizes L. Notice 
the logarithmic vertical scale. The vertical dashed line repre- 
sents Tc ^ 2.269. 



FIG. 5: (Color online) Acceptance probabilities A{T) as a 
function of temperature T for the ferromagnetic Ising model. 
The inset shows the acceptance rates as a function of temper- 
ature in the optimized case for varying system sizes L and a 
fixed number of temperature points. The vertical dashed line 
marks the critical temperature. 




FIG. 4: (Color online) Temperature sets for the ferromag- 
netic Ising model for different feedback steps. Starting from 
a geometric progression temperature set (step 0), we apply 
a feedback loop until the temperature set converges. While 
the geometric progression places many temperatures at low 
temperatures, the density of temperatures after the feedback 
optimization is highest at the bottleneck of the simulation 
around the critical temperature (marked by a vertical dashed 
line). Rapid convergence of the optimized temperature set is 
found after 3-4 feedback steps and a total of A'sw ~ 1.6 ■ 10^ 
swap moves in our parallel tempering simulations. 



replicas in the critical region by reallocating resources to- 
wards this temperature range. In contrast, the density of 
temperatures close to the lowest temperature is greatly 
reduced, thereby suppressing constant swapping of repli- 
cas at low temperatures where for the initial tempera- 
ture set multiple replicas converged to the fully-polarized 
ground state configuration. 

This effect becomes even more evident when measuring 
the acceptance probabilities for swap moves as illustrated 
in Fig.|Sl The acceptance probabilities for the initial tem- 
perature set based on a geometric progression saturate 
close to unity for temperatures below T < 0.7, whereas 
they show a pronounced dip already for small system 
sizes (L = 20) around the critical temperature (marked 
by a vertical dashed line). In contrast, the feedback- 
optimized temperature set shows a pronounced peak in 
the acceptance rate A(T) near the critical temperature 
where temperature points are accumulated by the feed- 
back. Away from the critical temperature region the ac- 
ceptance probabilities drop. The inset of Fig.[51shows the 
acceptance probabilities A(T) for the optimized tempera- 
ture sets for a fixed number of replicas and varying sizes of 
the simulated system. While the acceptance probability 
around the critical temperature remains nearly constant, 
the exact value away from the critical region decreases 
with increasing system size. This "mean" acceptance 
probability away from the bottleneck of the simulation 
can be tuned by varying the number of simulated repli- 
cas. 

The fact that for the optimized temperature set the 
acceptance probabilities vary with temperature contra- 
dicts various alternative approaches in the literature 
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FIG. 6: (Color online) Average round-trip times Trt before the 
optimization divided by the average round-trip times after the 
feedback optimization (t^^ *) as a function of system size. The 
data for the filled squares are for a system starting from a geo- 
metric progression and represent the speedup obtained by the 
feedback method. The open symbols correspond to a temper- 
ature set which initially has "flat" acceptance probabilities. 
The dashed lines are guides to the eye. 

[m IT^ IrR IT5 IT^ that aim at choosing a tempera- 
ture set in such a way that the acceptance probability of 
attempted swaps becomes independent of temperature. 
Naively, one might assume that the choice of constant 
acceptance rates produces an unbiased random walk in 
temperature space. This assumption is similar to the idea 
underlying generalized-ensemble algorithms that aim to 
sample a flat histogram in the energy such as the multi- 
canonical method 13, 0] or the Wang-Landau algorithm 
For these flat-histogram algorithms, it has been 
shown that they cannot reproduce the scaling behavior 
of an unbiased Markovian random walk in energy space, 
but experience critical slowing down 0. This slowdown 
can be overcome by optimizing the simulated statistical 
ensemble by a similar feedback algorithm ^ as presented 
here and sampling an optimized histogram that in general 
is not flat, but reallocates resources towards the bottle- 
neck of the simulation, e.g., in the vicinity of a phase 
transition or close to the ground state of the system. 

Measuring the diffusion of replicas in a subsequent sim- 
ulation for the optimized temperature set, we find that 
the current of replicas wandering from the lowest to the 
highest temperature is now characterized by a constantly 
decreasing fraction f{T) in agreement with Eq. ^ as 
shown in Fig. |21(triangles). In addition, we find that for 
the optimized temperature set the replicas wander evenly 
up and down in temperature space, that is 

'7~up ~ ^down ■ 

In Fig. El we show the ratio between the mean round- 
trip times Trt before optimization from a geometric and 
"flat" temperature set divided by the mean round-trip 



times T°f* after optimization in order to illustrate the 
speedup in replica diffusion attained by the feedback 
procedure. The data show clearly for all system sizes 
that the round-trip times after the optimization of the 
temperature set do not increase as fast as for a geomet- 
ric progression or "flat" temperature set. Note that for 
these temperature sets with a fixed number of tempera- 
ture points the average round-trip times before and after 
the feedback optimization scale ~ exp(aL^) for the sys- 
tem sizes studied. It is important to note that our algo- 
rithm identifies the bottleneck of the parallel tempering 
simulation that in this model occurs in the form of criti- 
cal slowing down at the phase transition solely based on 
measurements of the local diffusivity. The feedback then 
allows to shift additional resources towards this bottle- 
neck in a quantitative way. In the next section, we apply 
the algorithm to a more complex model with frustration 
and a different type of phase transition at zero tempera- 
ture. 



B. Fully- frustrated Ising model 

The Hamiltonian of the fully-frustrated (FF) Ising 
model is given by 

Hff = — ^ JijSiSj , (13) 

where the spins lie on the vertices of a two-dimensional 
square lattice with periodic boundary conditions. The 
bonds Jij are chosen such that \Jij\ = 1, but with the 
constraint that the product of the bonds around all pla- 
quettes of the system is negative, i.e., 

njy=-i. (14) 

□ 

The model does not order at finite temperatures, but ex- 
hibits a critical point at zero temperature. In the vicin- 
ity of this transition to a highly degenerate ground-state 
manifold, the system relaxes very slowly. 

In this section, we discuss our feedback optimiza- 
tion algorithm for two choices of the initial tempera- 
ture set. First we start from the temperature set in- 
troduced in Sec. IIV Al computed with a geometric pro- 
gression, Eq. 12)), which has T^in = 0.1, Tmax = 10.0 and 
M = 21 temperatures. In this first approach, we keep 
the number of temperature points constant for all sys- 
tem sizes L. In the second approach, we choose an initial 
temperature set where we fix again Tmin and Tmax to the 
aforementioned values but tune the number of tempera- 
tures M as well as their position in such a way that we 
obtain acceptance probabilities for swap moves that are 
approximately independent of the temperature ("flat") 
with a mean value of A{T) ~ 40% and deviations around 
this mean value of maximally 10% [2^ . 

We show the fraction /(T) of replicas diffusing from 
the lowest to the highest temperature as a function of 
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FIG. 7: (Color online) Fraction f(T) of replicas diffusing from 
the lowest to the highest temperature for the fully frustrated 
Ising model. Displayed are data for an initial "flat" tem- 
perature set with M = 21 temperature points that yields 
temperature-independent acceptance probabilities for swap 
moves (open squares). In addition, data for a geometric pro- 
gression (M = 21) are also shown (filled squares). After the 
optimization, the change in the fraction is independent of the 
temperature index (triangles). The inset shows the fractions 
as a function of temperature T. Data for A^sw = 2 ■ 10^. 



the temperature index in Fig. \7\ Similar to the Ising 
model, the data for the geometric progression tempera- 
ture set deviate considerably from a straight line which 
is expected for the optimal temperature distribution. A 
similar behavior is found when starting from a temper- 
ature set that initially has temperature-independent ac- 
ceptance probabilities. The local diffusivity in tempera- 
ture space calculated from the measured diffusive current 
is plotted in Fig. |S1 There is a pronounced dip in the dif- 
fusivity around T « 0.5 that we can identify as the tem- 
perature region where the system enters the highly de- 
generate ground-state manifold, e.g., by calculating the 
energy of the system, as plotted in Fig. ^ Again we 
find that the diffusivity points to a strong bottleneck of 
the simulation at the critical point which for the fully 
frustrated Ising model is at the transition to the zero- 
temperature ground-state manifold. The general shape 
of the diffusivity in the vicinity of this bottleneck re- 
mains unchanged with respect to the feedback. By apply- 
ing the feedback method, additional temperature points 
are shifted towards this bottleneck which is illustrated 
in Fig.^lfor the geometric progression (full symbols) as 
well as the "flat" temperature set (open symbols). For 
moderately large system sizes, we again find rapid con- 
vergence of the generated temperature sets within 3-4 
feedback steps and a total of Ngw ~ 7.5 • 10^ swap moves. 
For the optimized temperature set, the diffusive current 
is again characterized by a fraction of replicas drifting 



FIG. 8: (Color online) Local diffusivity D(T) of a random 
walk in temperature space for the fully-frustrated Ising model 
as a function of temperature T after the feedback optimization 
for several system sizes L. Notice the logarithmic vertical 
scale. 
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FIG. 9: (Color online) Energy per spin e{T) = {l/N)[{n)U 
as a function of temperature T for the fully-frustrated Ising 
model for several system sizes. The data show that already 
for T < 0.5 the energy is independent of temperature, thus 
signaling that the system has reached the ground state. The 
inset zooms into the temperature range around T — 0. 



from the lowest to the highest temperature that linearly 
decreases with the temperature index, see Fig. (trian- 
gles). 

The acceptance probabilities A{T) for replica swap 
moves are shown in Fig. 1111 for simulations with a ge- 
ometric progression and the optimized temperature set. 
While the acceptance probabilities peak at unity close 
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FIG. 10: (Color online) Temperature sets for the fully- 
frustrated Ising model for different feedback steps. Starting 
from a temperature set where the acceptance probabilities 
are independent of temperature with M = 21 temperature 
points (open symbols) and a temperature set obtained by a 
geometric progression also with M = 21 temperature points 
(filled symbols), we apply repeated feedback steps until the 
temperature sets converge to the optimal distributions. Also 
shown are data for an initial temperature set with M = 21 
equidistant temperature points (stars). Independent of the 
initial temperature set, an optimal temperature distribution 
is found after 3-4 iterations and ~ 7.5 ■ 10^ swaps. After 
the successful feedback, temperature points accumulate near 
the transition to the ground state slightly above zero temper- 
ature. 



to zero and dramatically drop in the geometric progres- 
sion temperature set, in the optimized set most tem- 
peratures are reshuffled to the low-temperature region 
slightly above zero temperature where the system enters 
the highly degenerate ground-state manifold. The inset 
to Fig. 111! shows the acceptance probabilities as a func- 
tion of temperature for a fixed number of temperatures, 
as well as T,„in — 0.1 and Tmax = 10.0 fixed. As in the 
case for the ferromagnet, the "mean" acceptance prob- 
ability away from the ground state bottleneck seems al- 
most independent of temperature and settles at a value 
that is determined by the number of temperatures used 
for a given system size L. 

In order to test the efficiency of the feedback method 
on the FFIM, in Fig. E|we show the ratio between the 
mean round-trip times Tj-t before optimization divided by 
the mean round-trip times r°^^* after optimization in or- 
der to illustrate the speedup in replica diffusion attained 
by the feedback procedure. The data show clearly for all 
system sizes that the round-trip times after the optimiza- 
tion of the temperature set do not increase as fast as for 
a geometric progression or "flat" temperature set. For 
these temperature sets where the number of temperature 
points increases with system size we find that the average 



FIG. 11: (Color online) Acceptance probabilities A{T) for 
replica swap moves as a function of temperature T for the 
fully-frustrated Ising model. While the acceptance probabili- 
ties for a geometric progression temperature set show a pro- 
nounced dip close to T = 0, the optimized ensemble shows a 
peak close to zero temperature where the system enters the 
ground-state manifold. The inset shows, for a fixed number 
of temperatures, the acceptance rates as a function of tem- 
perature for different system sizes L. As for the Ising model, 
the "mean" value away from the bottlenecks can be tuned by 
increasing the number of temperatures. This illustrates that 
in order to obtain higher acceptance rates away from the bot- 
tlenecks of the simulation, the number of temperatures have 
to be increased with increasing L. 

round-trip times scale a + bx'^. 

Finally, we discuss the effects of varying the number of 
temperatures M in the temperature set. Figure lOl shows 
the average round-trip times for the fully frustrated Ising 
model (L = 20) as a function of the number of tempera- 
tures M . For M > 12, the average round-trip times show 
only moderate variations with the number of temperature 
points M, whereas for a smaller number of temperatures 
the average round-trip times increase drastically. This 
can be understood by keeping in mind that a parallel 
tempering swap will only be accepted with a high prob- 
ability if the energy distributions between neighboring 
temperatures overlap. If the minimum and maximum 
temperatures are fixed and M is reduced, the energy dis- 
tributions will cease to overlap, which accounts for the 
increased average round-trip times. Factoring in the total 
CPU time, which we define as the product of the aver- 
age round-trip time with the number of temperatures, 
the minimum is more pronounced (inset to Fig. I13|l and 
clearly illustrates that while a larger number of temper- 
atures has little effect on the round-trip times, the total 
CPU time increases drastically with increasing M . 

Because the data for the average round-trip times vs 
M have an optimal value, it is conceivable to develop a 
feedback optimization method that optimizes both the 
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FIG. 12: (Color online) Average round-trip times Trt before 
the optimization divided by the average round-trip times after 
the feedback optimization (t^^'*) as a function of system size 
for the fully-frustrated Ising model. The data for the filled 
squares are for a system starting from a geometric progres- 
sion temperature set and represent the speedup obtained by 
the feedback method. In addition, we show data for a tem- 
perature set with "flat" temperature-independent acceptance 
rates (open squares). The dashed lines are guides to the eye. 




FIG. 13: (Color online) Average round-trip times T^t as a 
function of the number of temperatures M for the fuUy- 
frustrated Ising model with L = 20 after the feedback op- 
timization. The data show that the round-trip times only de- 
pend moderately on the number of temperatures M, provided 
that there is sufficient overlap of the energy distributions. For 
a small number of replicas, this is no longer the case and the 
round-trip times increase drastically, e.g., for M < 12 in this 
plot. The inset shows the CPU time which is the product of 
the average round-trip time with the number of temperatures 
M. The data show a more pronounced minimum. 



position of the temperatures and the number of temper- 
atures M. Furthermore, in addition to optimizing the 
number and locations of the temperature points, we have 
also explored varying the ratio of the number of sweep 
moves attempted to the number of replica-swap moves 
attempted, since this is yet another parameter one must 
set in a parallel tempering simulation. This ratio can be 
adjusted globally (the same ratio at all temperatures) or 
locally (the ratio optimized independently at each tem- 
perature). This will be discussed in more detail in a 
subsequent communication. 



C. Spin glasses 

Since the optimization of temperature sets improves 
the sampling for the two paradigmatic spin models dis- 
cussed above, it is a natural step to ask how this feed- 
back optimization technique can be applied to improve 
statc-of-thc-art parallel tempering simulations of Ising 
spin glass models, such as the three-dimensional (3D) 
Edwards- Anderson Ising spin glass ,21,, j2§| : 



I J Si Sj . 



(15) 



Here the spins lie on the vertices of a cubic lattice with 
periodic boundary conditions. The bonds Jij are chosen 
according to a Gaussian distribution with zero mean and 
standard deviation unity. The system undergoes a spin- 
glass transition at = 0.951(9) HHHm. 

For the spin glass there is the additional complexity 
that different disorder realizations can lead to strong 
sample-to-sample variations. Thus one could also sur- 
mise that strong sample-to-sample variations exist for 
the time it takes to equilibrate individual samples. This 
becomes evident when measuring the round-trip times 
for replicas in a standard parallel tempering simulation 
with a fixed temperature set, as illustrated in Fig. [HI for 
the Edwards- Anderson Ising spin glass. We find that the 
measured round-trip times follow a fat-tailed Frechet dis- 
tribution [sol l3ll| (solid line, fit performed with R 32j). 
The integrated generalized extreme value distribution is 
given by: 



H^-f,-p{T) = exp 



- 1 + 



(16) 



Here fj, represents a generalized most probable value (lo- 
cation parameter) and /3 a standard deviation (scale pa- 
rameter) . The value of the shape parameter ^ determines 
whether the distribution is thin-tailed (f < 0, WeibuU), 
Gumbel (C = 0), or fat-tailed (^ > 0, Frechet) ^QJ. Note 
that when f > 0, the m-th moment of the Frechet dis- 
tribution exists only if |^| < 1/m, e.g., if ^ > 1/2 the 
variance of the distribution is not properly defined. Our 
results are in agreement to similar observations for broad- 
histogram simulations 0, IsM Is^ . Note that the distri- 
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FIG. 14: (Color online) Distribution of average round-trip 
times for 5000 different samples of the 3D Edwards- Anderson 
Ising spin glass with Gaussian disorder and fixed system size 
L = 4 in the temperature range from Tmin = 0.10 to Tmax = 
2.0. The data follow a fat-tailed Frechet distribution (solid 
line) with a shape parameter ^ — 035(1). The inset shows the 
shape parameter ^ as a function of system size L. Already 
for L > 5, the shape parameter becomes ^ > 1/2, indicating 
that the variance of the distribution is no longer properly 
defined. The simulations have been performed using a fixed 
temperature set with M — 27 temperature points distributed 
such that, on average, replica swap moves have a nearly flat 
acceptance rate. 



bution is increasingly more fat-tailed as the system size 
increases (see the inset to Fig. I14|l . 

The measurement of the round-trip times thus allows 
to classify individual spin-glass samples as "typical" with 
round-trip times in the bulk of the distribution or "hard" 
with round-trip times in the tail of the distribution. Such 
a classification can be very useful to shift computational 
resources towards the "hard" samples in the course of 
a simulation as these samples potentially might require 
substantially longer simulation times in order to equi- 
librate. Preliminary tests suggest correlations between 
round-trip and equilibration times. The observation of 
strong sample-to-sample variations in the distribution of 
round-trip times also raises the general question, whether 
previous spin-glass studies have properly equilibrated the 
"hardest" samples in their simulations. It remains to be 
verified whether this introduces a systematic error in the 
analysis of spin-glass systems. Specifically, the finite-size 
scaling should be sensitive to such systematic errors, as 
it has been observed that the number of "hard" samples 
significantly increases with system size, see the inset in 
Fig. HI and Refs. 2, 33, andlM 

On the other hand, we can ask whether we can opti- 
mize the simulated temperature set and generate a "com- 
mon" temperature set for samples from the various parts 



of the distribution. To do so, we apply the feedback op- 
timization outlined above in such a way that we generate 
a common probability distribution fjiT) for a set of sam- 
ples that are each characterized by their own diffusivity 
Di{T), steady-state fraction fi{T) and current ji, which 
are related by 

j,=D,iT)fjiT)^, (17) 

where the index i indicates the samples in the given set. 
Rearranging this equation, the local diffusivity of an in- 
dividual sample can be expressed as 

To ensure equilibration of all samples we want to simulate 
each sample for a fixed number of round trips, despite the 
strong sample-to-sample variations. In order to minimize 
the overall computer time spent to simulate such a set of 
samples, we then minimize the sum of round-trip times 
^^Ti. This is equivalent to minimizing the sum of the 
inverse of all currents, i.e., , since the current ji is in- 
versely proportional to the round-trip time . Following 
a similar line of arguments as for the original tempera- 
ture set optimization, we find that the optimal common 
temperature distribution 77°p*(T) is proportional to the 
square root of the sum of inverse local diffusivities 

Again we can use a feedback loop to find an optimized 
common temperature set by feeding back the measured 
local diffusivities 

fi°''\T)^C^fj{T)J2{r,-^^ , (20) 

where C is a normalization constant. The common opti- 
mized temperature set is then found using a partial inte- 
gration as given in Eq. Hll|l. 

V. CONCLUSIONS 

We have introduced an approach to systematically op- 
timize temperature sets for parallel tempering Monte 
Carlo simulations using an adaptive feedback method 
that is motivated by a recently developed ensemble op- 
timization technique for broad-histogram Monte Carlo 
simulations 0. We have applied the method to two 
paradigmatic spin models, the ferromagnetic Ising model 
and the fully frustrated Ising model in two dimensions. 
For both models we have shown that the feedback tech- 
nique improves sampling of the phase space by reducing 
the average round-trip time of replicas diffusing in tem- 
perature space. 
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Probably our most important result is the insight, 
that the common wisdom that temperature sets in par- 
allel tempering Monte Carlo should yield temperature- 
independent acceptance probabilities for the swap moves 
is not necessarily an optimal choice. Our feedback al- 
gorithm shifts temperature points in the optimized tem- 
perature sets towards the bottlenecks of the simulation, 
typically in the vicinity of a phase transition, where equi- 
libration is suppressed. In particular, this has the effect 
that the acceptance probabilities for replica swap moves 
are higher around the bottlenecks and are not tempera- 
ture independent for the so-optimized temperature sets. 

We also outline an approach to define "common" tem- 
perature sets for systems that require configurational 
averages, such as spin glasses. In addition, we have 
briefly discussed the effects of sample-to-sample fluctu- 
ations with respect to equilibration times of individual 



spin glass samples, thus pointing towards a potential 
source of systematic errors in previous numerical stud- 
ies of spin glasses. 

Clearly a deeper analysis of feedback-optimized paral- 
lel tempering Monte Carlo is needed, in particular in the 
context of spin glasses, as well as the questions raised at 
the end of Sec. HVBl 
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